Analytical modeling for the heat transfer in sheared flows of nanofluids 
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We developed a model for the enhancement of the heat flux by spherical and elongated nano- 
particles in sheared laminar flows of nano-fluids. Besides the heat flux carried by the nanoparticles 
the model accounts for the contribution of their rotation to the heat flux inside and outside the 
particles. The rotation of the nanoparticles has a twofold effect, it induces a fluid advection around 
the particle and it strongly influences the statistical distribution of particle orientations. These 
dynamical effects, which were not included in existing thermal models, are responsible for changing 
the thermal properties of flowing fluids as compared to quiescent fluids. The proposed model is 
strongly supported by extensive numerical simulations, demonstrating a potential increase of the 
heat flux far beyond the Maxwell-Garnet limit for the spherical nanoparticles. The road ahead 
which should lead towards robust predictive models of heat flux enhancement is discussed. 



I. INTRODUCTION 

The enhancement of heat transfer by embedded nano- 
particles in a fluid subjected to a temperature gradient 
is an important issue that is expected to help in techno- 
logical application like solar heating, and various cooling 
devices, including miniaturized computer processors. Ac- 
cordingly, many experiments were conducted in the last 
couple of decades [iJ-Q , with a breakthrough announce- 
ment from a group at Argonne National Laboratory, who 
studied water and oil-based nanofluids containing cop- 
per oxide nanopaticles, and found an amazing 60% en- 
hancement in thermal conductivity for only a 5% volume 
fraction of nanoparticles @. Subsequent research has 
however generated what Ref. referred to as "an as- 
tonishing spectrum of results" . Results in the literature 
show sometime enhancement in the thermal conductiv- 
ity compared to the prediction of the Maxwell-Garnett 
effective medium theory, and sometime values that are 
less than the same prediction [ll-dl ■ Confusingly enough, 
these discrepancies occur even for the same fluid and the 
same size and composition of the nano-particles. 

Some of these conflicting results were explained by ei- 
ther the formation of percolated clusters of particles (for 
the case of enhancement with respect to the Maxwell- 
Garret prediction) or by surface (Kapitza) resistance (for 
the case of reduction with respect to the same prediction) 
Q. Interestingly enough, enhancement is typically seen 
in quiescent fluids, and one expects that the agglomera- 
tion of clusters will become impossible in flows, be them 
laminar or turbulent. Of course, in technological appli- 
cation flowing nano-fluids may be the rule rather than 
the exception. Thus our aim in this article is to develop 
models of heat flux in flowing nano-fluids, where the flow 
can be either laminar or turbulent. In such systems we 
cannot expect an enhancement of heat flux due to the 
agglomeration of particles (at least at low volume frac- 



tions), and for the sake of simplicity we will assume that 
there is no Kapitza resistance. Since it was reported in 
the literature that elongated nano-particles are superior 
to spheres in quiescent nanofluids |8| , we will study both 
spheres and spheroids in flowing nanofluids. We will ar- 
gue that gencrically elongation may not be advantageous 
at all, and will explain why. 

For completeness, we will begin by studying quiescent 
nanofluids under a temperature gradient. We will present 
extensive numerical simulation that will demonstrate a 
very good agreement with the Maxwell-Garnett theory 
for spherical nano-particles and with the generalization 
of Nan et al for elongated particles. In the case of 
flowing nanofluid we will offer a model that will provide 
expressions for the heat transfer for different values of the 
aspect ratio of the particles, for different volume fractions 
and for laminar and turbulent flows. 

The structure of the paper is as follows. In the next 
Section (|TT| we formulate the problem, describe the equa- 
tions of motion and the numerical procedure, and discuss 
the dilute suspension approximation. In Section IIIII we 
describe the results of numerical simulations of spherical 
and spheroidal particles in a fluid at rest and in a shear 
flow. In Section IIVI we develop an analytical model of 
heat transfer characteristic, i.e. the Nusselt number, and 
analyze the model predictions of heat flux enhancement 
in two limiting cases of very strong Brownian diffusion 
and a weak one. 



^ A spheroid, or ellipsoid of revolution is a quadric surface obtained 
by rotating an ellipse about one of its principal axes; in other 
words, an ellipsoid with two equal semi-diameters. 
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II. NANOFLUIDS IN THE DILUTE LIMIT 

In this section we formulate the model for dilute 
nanofluids laden with elongated spheroidal nanoparticles. 
This includes the basic equation of motion for the ve- 
locity and temperature fields and the boundary condi- 
tions (with constant velocity and temperature gradients 
far away from particles). We also present details of the 
numerical simulations and their validation. 



A. Formulation of the problem and the flow 
geometry 

When the nanofluid is very dilute one can disregard 
the effect of one particle on the other and consider the 
nanofluid as an ensemble of noninteracting particles. One 
of these is shown in Fig. [TJ The center of the spheroidal 
particle is at the position x~y = z = Om the middle 
of a plane Couette flow between two _ff-separated hor- 
izontal parallel walls that move in the i-direction with 
opposite velocities V± = ztV/2. Due to the symmetry, 
the forces are balanced, and the particle neither migrates 
nor collides with the walls; the particle's center remains 
at (0,0,0). This allows us to eliminate in the numerics 
any effect of the particle sweeping parallel to the walls. 
Note that this effect is anyway absent in homogeneous 
cases. 

We employ periodic boundary conditions along the x 
and the z direction. This results in a periodic replication 
of the computational box (together with the particle) in 
the XZ -plane, leading to a "monolayer" of an infinite 
number of periodically distributed particles in the XZ- 
plane. Remarkably enough (and by reasons that will be 
explained below), this allows us to reproduce in numerics 
with a single particle effects of a finite volume fraction tp 
at least up to iy9 < 20 30%. 

The walls are kept at fixed temperatures T± = T±A/2. 
In the absence of particles these boundary conditions give 
rise to a vertical temperature gradient V^T = A/ H and 
shear S = S^y = V/H (see Fig. d]). 

The spheroidal particle's semi-axes are a (the longest) 
and b (the shortest). a>b (Fig. [1]). The thermal diffu- 
sivity inside the particle is Xp. The carrier fluid has a 
kinematic viscosity v and a thermal diffusivity Xf- For 
simplicity, at this stage we neglect the effect of gravity 
and the particle mass. The co-ordinate system and polar 
angles are given by: 

rix = cos 6 =sin0sin0, (la) 
Tiy = sin 6 sin (p = sin 6 cos (j> = cos 9^ , (lb) 
71-2 = sin cos = cos , (Ic) 

where rii is the projection of the unit- vector n onto the 
axis i, i = {x,y,z}, and is the angle between the 
particle's largest axis and the j/-axis. 

In the absence of particles, changing the intensity of 
the laminar shear fiow does not change the heat fiux since 




FIG. 1. The particle in a plane Couette flow. The largest 
particle axis is shown as a (red) rod. The velocity field (far 
from the particle) is {Sy, 0, 0) where the shear S = (V+ — 
V-)/H in the particle-free flow with H being the distance 
between the walls. For the theoretical development we assume 
the creeping flow conditions as in [9l-[lll|. 



there is no velocity orthogonal to the wall and the system 
is homogeneous in the direction parallel to the wall. In 
the presence of a particle, the shear induces it to rotate, 
the simplest case being a spherical particle which rotates 
at a constant angular velocity 0,^ — S/2 (provided the 
particle radius is much smaller than the distance to the 
wall). This rotation induces a vertical fiow motion near 
the particle. The modification of the velocity profile has 
three consequences. The first directly affects the heat 
convection in the fiuid through a convective contribution 
VyT. The second comes from the particle rotation which 
brings up the hotter side of the particle during its rota- 
tion. The third changes the heat conduction due to a 
modification of the temperature profiles at the surface of 
the particle. 



In the case of non spherical particles the presence of a 
shear has also a strong infiuence of the statistical distri- 
bution of particle orientations. 
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B. Dimensionless parameters 

To fix the notation wc fist here the most important 
dimensionless parameters involved in the problem: 



aspect ratio: 


r 


= a/6, r>l, 


(2a) 


diffusivities ratio: 


k 


= Xp/Xf , 


(2b) 


Prandtl nmnber: 


Pr 


= v/xt , 


(2c) 


Reynolds number: 


Re 


= {2afS/v, 


(2d) 


Peclet number: 


Pc 


= {2afS/xv ■ 


(2c) 


Note that the particle's lar{ 


^est 


axis defines our 


length- 



scale. Also, 



Pe = RePr/fc . 



(3) 



C. Basic equations of motion 



The dynamical effects previously mentioned can be 
quantitatively studied by solving the following system of 
equations for the temperature in the fiuid, Tf , and in the 
particle, Tp, 

^ + («.V)Tf = XfV2rf, ^ = XpV%. (4a) 

together with the boundary conditions, 

at the surface: Tf = Tp , XfVTf = XpVTf , 

on the walls: Tf (a;, 0, 2) = T„ , Tf (cc, H, z) = T+ . 

The fluid velocity in the whole domain can be found 
by solving the incompressible Navier-Stokes equations: 



dt Pi 



vV^Uf, V-Uf = 0, (4b) 



together with non-slip boundary conditions at the surface 
of the particle and at the walls: 

Mf(x,0,z) = {V-,0,0}, Uf{x,H,z) = {l/+,0,0}, 

ttf (x, y, z) = Mp(a;, y, z) at the particle's surface. 

Here pf is the carrier fiuid density and p is the pressure. 
Clearly, on nano-scales the temperature difference (across 
a particle) is sufficiently small to allow us to neglect the 
dependence of the fluid's and particle's material param- 
eters on the temperature, i.e. we take Xf ^^id Xp 
constants. 

The dynamics of a neutrally buoyant particle isgov- 
erned by the equations of the solid body rotation [T^]: 





^ X 






h 


dt 






Ice 




dd^^ 










dt 










d^t^ 


T?) 









dt 



^^(b)^(b) 



(4c) 



where Cl^^'' is the angular velocity in the body- fixed 
frame, T^'"^ is the torque in the body-fixed frame and 
/ is the moment of inertia tensor. 
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FIG. 2. Temperature profile along a line passing through 
the center of a spherical particle (of radius R) in a quiescent 
fluid. Dashed (green) line - conductive temperature profile 
in the pure fiuid, solid (red) curve - Eq. ([ASP jlil ] for an 
infinite domain with a temperature gradient, dots - numerical 
results (at lattice points) with k = 5 (full circles) and k = 100 
(squares). The present test was performed with Pr = 1, H = 
64 and i? = 10 lattice units. 



D. Numerical approach and its validation 

The numerical simulation of the conjugated heat trans- 
fer problem given by Eqs (|4]) is performed by means of 
two coupled D3Q19 Lattice Boltzmann (LB) equations 
under the so-called BGK approximation Details of 
the simulations are presented in Appendix |X] . Besides 
purely numerical means, the control of the simulations 
was done by its comparison with known analytical solu- 
tions. As an example. Fig [5] shows the comparison be- 
tween the analytical temperature profile (|A3p (solid red 
line) and numerical profile (black dots) for a moderately 
big spherical particle [H/{2R) = 3.2]. The second test 
was the comparison in Fig. [TU] of the particle's angular 
velocity theoretically obtained by Jeffery Q and simu- 
lated by our LB-realization. For more details, again see 
Appendix [X] 



E. From single particle modeling to finite volume 
fraction 



Here we discuss how to relate the contribution of a sin- 
gle particle to the heat flux with the total contribution of 
many weakly-interacting particles randomly distributed 
in the flow occupying a finite volume fraction Lp. 

The first step is to consider fully dilute limit, </? — > 0, 
in which we can exploit the fact that the particle contri- 
bution is additive and proportional to ip. In this way we 
can introduce the Nusselt number as the ratio between 
the total heat fiux (inside the computational cell) and 
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the conductive heat flux in the basic ceU: 



NU : 



cond/ 



-XfVyT 



XeS 

Xf 



(5) 



where q is the total heat flux, (/cond = "Xf'^yT is the 
conductive part and Xeff is the effective heat difFusiv- 
ity of the composite fluid. Also, V^T = dT{x,y, z)/dy. 
Without particles (q) = (gcond) and Nu = 1. For dilute 
suspensions, in which ip <ti 1, Nu can be expanded in 
powers of ip p^ : 



Nu 



Ki + K2(p- 



(6a) 
(6b) 



where Ki are dimensionless constants dependent on the 
Peclet number, on the particle's aspect ratio, r, on the 
relative heat diffusivity, fc, and maybe other parameters 
like Reynolds or Prandtl numbers [sec Eqs. ©J, i.e. Kj ~ 
if,(Pe,r,fc,...).^ 

As expected, in our simulations the quantity (Nu— 1) 
is indeed proportional to ip for (p ^ 0.05. This allows us 
to find Ki as a function of other parameters ([2]) of the 
problem. 

The next order term in the expansion (j6bp . i.e. K2 4>, 
contains very important information about how Nu de- 
pends on ip for small but finite values of tp. Generally 
speaking, to get this information from numerics one needs 
to solve Eqs. (|4]) for many particles, randomly distributed 
in space. We demonstrate that this problem can be 
tremendously simplified by reducing it to a one-particle 
case, in which this particle of volume Vp = (pVceii is put 
in the center of computational box of volume Vcdi with 
periodic boundary conditions in the horizontal (orthog- 
onal to the temperature gradient) x- and z-directions. 
In this way the total system can be considered as con- 
structed from a periodic repetition of the basic elemen- 
tary cell in the x- and z-directions. Comparing in Sub- 
sec. [111^1] our numerical results with analytical findings 
obtained under the assumption of random particle dis- 
tributions, we see that the precise particle distribution is 
not essential - what really matters is the actual volume 
fraction. 

The reason for that is quite simple: as one sees in 
Fig. [21 the deviation from the linear temperature profile 
becomes important at distances^ smaller than the par- 
ticle diameter 2R for both fc = 5 and k = 100. Thus 
any nonlinear dependence of Nu on ip can appear only 
if the particles are sufficiently close to overlap the devi- 
ation from linear temperature profile. We see from our 
numerics that the this interaction causes 10% deviation 
from the straight line in Nu vs. tp dependence for vol- 
ume fractions of about 0.1 and more for k < 100. For 



^ It can be shown from Eq. IIA3I I 

smaller than R^5~^(k — l)/(fc -I- 2)|^''', where 5 is the deviation 
criterion, i.e. ITf/Tnn — 1\ < S and Tn^ = Gy 



that the distance should be 
11/3 
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FIG. 3. (Color online). Quiescent fluid, spherical particle in 
the center. Nu((p). Dashed curves - theory, Eq. (O and dots 

- simulations for different k (color-coded from "cold" -blue, 
fc = 1, to "hot"-red, fc = 100). Insert: Ki{k). Dashed curve 

- theory, Eq. ((Tjl and dots - simulations for different k (same 
color-code). Notice, that for 128"^ runs the accuracy is better 
then 0.3%, while for 64'^ runs it is about 0.7%, and only for 
really small 32^ runs it reaches the worst value of 1.5%. These 
means that for our purposes, usually the runs with 64^ box is 
fully satisfactory. 



ip > 0.1 there is some nonlinear dependence of Nu on ip 
but it practically coincides for the random and the peri- 
odic distributions of particles. 



III. TOWARDS AN ANALYTICAL MODEL OF 
THE HEAT FLUX 

In this Section we begin with collecting the required 
information about the dependence of Kj on the various 
parameters To do this we consider simple limiting 
cases, for which some of these parameters are put to zero. 
We start with the case of spherical particles. 



A. Spherical nano-particles 

1. Test case: s spherical nano-particle in a fluid at rest 

Here we consider the simplest case of a spherical par- 
ticle (r = 1) in a fluid at rest (Re = Pe = 0). Chiew and 
Glandt [lB| suggested the following formula for this case: 



Kik,p) 



3Ki{k) 



Nu = 1 



K,{k)^ ' 

3Ki{k)ip 
3-Ki{k)ip 



Ki{k) = 3 



fc- 1 
k + 2 



(7a) 
(7b) 



We notice that the expression for A'i(fc) which controls 
the very dilute limit {ip — ?► 0) is the same as in the 
Maxwell model [l3|. Obviously, for fc = 1 one has zero 
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FIG. 4. (Color online). A spherical particle in a shear (Couette) flow. Upper left panel: k — 1, ip = 9.3%, 14.8% and 19.4%. 
Upper right panel: k = 5 for wider region of ip, but Pe is limited by 10. Both upper panels show Nu(Pe)/Nuo — 1 vs. Pe, 
and both insets show Nu(Pe) vs. Pe. Here Nuo =Nu(0). Dashed (red) line is a fit oc Pe^, cf. Eq. (|8a|l : (color) symbols are 
simulations with different volume fractions: tp = 1.6%, 4.4%, 9.3%, 14.8%, 19.4%. Solid (black) curve on upper left panel 
represents the fit by Eq. (|9}. Lower panels: ip — 9.3% and different k (from "cold"-blue, fc = 1, to "hot"-red, k = 20). Dashed 
("temperature" -colored) lines are fits to simulated data (symbols, same colors) by Eq. (|8a|) . left panel, and by combined Eqs. (jSj 
and ([7]), right panel. Inset on lower left panel; B{k,ip) vs k at ip = 9.3%. Dot-dashed (brown) line is the fit to Eq. (jSbJ. 
Parameters: 64"^, Pr = 1, thus, Re — Pek. 



effect, i.e. Nu = 1. For fc oo, A'l = 3 and one has 
maximal possible enhancement (at fixed ip 1). For 
k = 2 one has 25% of this value, fc = 4 gives already 
50% of the effect, fc = 10 ^ 75%, fc = 20 ^ 86% and 
fc = 100 ^ 97%. Clearly, the larger k the better, but 
increasing k above 100 is not effective. 

Actually, Eq. ([7|) includes also information about the 
nonlinear dependence K{ip). However this dependence is 
very weak in the relevant range of parameters: e.g. for 
fc = 10 and iy9 = 0.2 the difference between K{ip) and 
its linear approximation is only 15%. For smaller k this 
difference is even smaller. The physical reason for that 
was explained at the end of previous Section. 

We tested the prediction of Eq. ([7]) by simulating the 
heat flux in a periodic box as explained above in Sec. |TT] 
and in Appendix For the fluid at rest the (volume 
averaged) Nusselt number (as a function of if at different 
values of k) was measured and is shown in Fig. [31 where 
the inset shows Ki{k). The overall conclusion is that 

Eq. ^ agrees extremely well with all our simulations and 



thus can he used in modeling the effect of spherical par- 
ticles in a fluid at rest. 



2. Spherical nano-particle m a shear flow 

The next important question is how the particle ro- 
tation affects the heat flux. As already mentioned, this 
rotation induces fluid motions around the particles thus 
causing convective contribution to the heat flux. This 
changes the temperature profile around the particles and, 
in turn, affects the heat flux inside the particles. To study 
these issues analytically is extremely difficult and thus 
numerical simulations can play a crucial role. Due to the 
fact that we have to deal with a number of parameters, 
we will carefully examine them starting with the sim- 
plest case of spherical particles. Results for spheroidal 
particles are discussed in Sec. IIIIB 21 

We begin by considering the case fc = 1; such a particle 
at rest docs not affect the heat flux, thus Nu — 1. In 
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Fig.m upper left panel, we present Nu(Pe)/Nu(0) — 1 as a 
function of Pe for k = 1 at various volume fractions (from 
9.3% to 19.4%), while the inset is showing Nu vs. Pe. 
Observing the data, we suggest for Pe^ 10 the following 
ansatz 

Nu(Pe) ~ Nu(0) [l + B{k, ip)Pe^] , (8a) 

shown in this panel as a straight dashed (red) line. One 
sees that Eq. ((8a|) is approximately valid up to Pe ~ 10, 
where the deviation is about 1%. We note also that for 
k = 1 the value of B is very small, B{l,ip) ~ 10~^. 
A similar analysis for k ~ 5, (cf. Fig. SI upper right 
panel), shows that Eq. (|8al) fits the data even better with 
a similar value -8(5, (p) ~ 2.5 • lO^"'. 

We see that B depends weakly on ip but more strongly 
on k. Thus, to determine the leading fc-dependence we 
put a sphere of radius 18 in a 64"^ computational volume 
(all in lattice units); this is equivalent to ip ^ 9.3%, see 
Fig. m lower panels. The fc-dependence of B{k, ip) in 
Eq. ([8a|) in the range of </? ~ 10% can be fitted by: 



S(fc,(p)~Bi((^) cxp[B2((p)fc], (8b) 
Bi ~ 9xl0"^ B2 ~ 0.15 . (8c) 

For larger Pe> 10 , the Nu vs. Pe dependence deviates 
down, as expected, because at Pe — >■ 00 this dependence 
should saturate. Our analysis shows that this dependence 
can be fitted with a good accuracy by the formula that 
generalizes Eq. ([8a|) : 



Nu(Pe) 
Nu(0) 



- 1 



B{k,(p) Pe^ 



1 + Pe/Pei(fc) + [Pc/Pe2(fc)] 



(9) 



The upper left panel in Fig. |3] shows by a solid (black) 
curve how this model works for fc = 1 (with Pei(l) w 61 
and Pe2(l) « 63). 

Having in mind that in many applications related to 
nanofluids the value of Pe is smaller than 0.01 and in 
any case rarely exceeds unity, we reach the conclusion 
that the convective heat flux around spherical particles 
and variations of the heat flux inside spherical particles 
due to their rotation can be neglected. Equation ([7|) can 
be used to model the heat flux enhanced by spherical 
nano-particlcs with finite volume fraction up to 25% and 
any actual value of k in fluids at rest and in shear flows. 

The next question to consider is the effect of the par- 
ticle shape. We will begin with the case of spheroidal 
particles in fluids at rest. 



B. Spheroidal nano-particles 

1. Spheroidal nano-particle in a fluid at rest 

The general expectation is that spheroidal nano- 
particlcs should be able to enhance the heat flux much 
more effectively then spherical particles. To achieve this 
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FIG. 5. (Color online) Spheroid in a quiescent fluid at dif- 
ferent angles Sn. Both panels: fc = 5, Pr = 1, 64"^. Symbols - 
simulations, solid (red or blue) curves - Eq. Ulip with appro- 
priate r and p, dashed (black) curves - fits by Eq. (|10p to the 
simulations. Upper panel: r — 2. Data: p ~ 10.79% (up- 
per), p ~ 4.25% (middle) and p ~ 0.69% (lower). The largest 
mismatch between the simulations and Eq. (jlip is 2.8% at 
p = 10.79%. Lower panel: lower (red) dots: p = 4.25%, 
r — 2.0, and upper (blue) squares: p = 4.31%, r = 2.7. 
The dotted (magenta) curve is Eq. (jlip with p — 4.25% and 
r = 2.7. The effect of changing r is larger than that of p. 



they have to be oriented in the "right" direction, i.e. 
along with the temperature gradient. Therefore, logi- 
cally, the study of the heat flux in nanofluids laden with 
spheroids should begin with the clariflcation of the ef- 
fect of the spheroid orientation on the heat flux. In this 
Section wc consider this effect for a given and fixed ori- 
entation of the spheroids. For this purpose we perform 
numerical simulations of the heat flux with spheroidal 
nano-particles in a quiescent fluid (Pe= 0) with a temper- 
ature gradient, see Fig. [U taking for concreteness fc = 5. 

For different orientations of the spheroid (i.e. different 
~ the angle between the temperature gradient vec- 
tor and the largest spheroid's axis), we measure different 
Nusselt numbers: when 9^ = spheroids with higher 
conductivities tend to forms a thermal shortcut, thus Nu 
is larger than for ~ 7r/2, where this shortcut effect is 
reduced (Fig. [5]). Our numerics shows that the dcpen- 
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dence of the Nusselt number on the spheroid orientation, 
i.e. Nu(0n)j for small (p is well described by a simple 
formula 



NU: 



1 - 



in ~^ ^aniD COS (6 



(lOa) 
(10b) 



and both /Cmin and /Camp are functions of fc, r and ip. 
This equation is a generalization of Eg. ^ for the case 
of spheroids, and it was suggested in [18f for the limiting 
case of r — > cx) (nanotubes) . For large (p > 5% deviations 
become visible. For finite if there have been attempts 
to study the role of the particles shape and form- factors, 
e.g. by Nan et al. [1]: 



IC 



/3i + (/32-/3i)(cos20n) 



1 - p [Li/3i + {L2I32 - Li/3i)(cos20N>] 
where: = [L^ + l/(fc - 1)]"^ 
L2 = 1 - 2Li , 



(11) 



r > 1: Li 



2(r2 - 1) 



arccosh(r) 



1 



Here (...) stands for an ensemble average. Clearly, in 
the case of a single particle as in our simulations, or a 
periodic array of the particles, the average coincides with 
the single value. Notice, also that for spheres, i.e. when 
r -> 1, Li = L2 = 1/3, /3i = /32 = Ki from Eq. (0), 
and Eq. pT|) simplifies to Eq. ([7]). For small (p Eq. PT|) 
coincides with Eq. (jlObp . giving analytical expression for 

/Cniin and /Camp- 

The range of validity of equation dTTI) was tested by 
means of a series of simulations at varying angles and 
volume fractions. The results are reported in Figure [S] 
We can conclude that the model in Eqs. PT|) well repre- 
sents the heat flux in the presence of resting spheroidal 
particles in the relevant range of the parameters k, r and 
p. To apply Eqs. ([TT|) for spheroids in a simple shear 
flow wc have to find how the ensemble average (cos^ 0n) 
depends on r. For this purpose we first need to know 
the orientational distribution function of spheroids in the 
shear fiow. This is a subject of the following Subsection. 



2. Spheroidal nano-particle in a shear flow 

The effect of rotation of elongated spheroidal particles 
(r > 5) is very similar to that of spherical ones (r = 
1), just the parameters i?, Pci and Pe2 of the advanced 
fit ^ depend on the aspect ratio r. As an example, we 
presented in Fig. [5] a preliminary result of the computed 
and fitted Nu(Pe) dependence for r = 5 and k — 5. In 
this case B fa 5.3 x 10-^ Pei « 5.2, Pe2 ~ 30.3. This 
is interesting to compare with the respective parameters 
for r = 1 and fc = 5: B « 2.5 x 10"'*, Pci « 61, Pci « 
63. One sees that the Pe-enhancement of the heat fiux 
with elongated particles is smaller then for spherical ones 
[B{r > 1) < B{r = 1)] and it saturates at smaller Pc: 
Pei(r > 1) <Pci(r = 1) and Pe2(r > 1) < Pe2(r = 1). 



Moreover, the overall conclusion that Nu(Pe) is almost 
Pe-independent for Pe^ 1 remains valid, thus, the model 
developed above in Eq. pT|) with Eq. ([^ may be safely 
used for nano-particles laden flows at Pe< 1. 




2 X 10"' 



2 3 5 10 

Peclet number, Pe 



FIG. 6. (Color online). Example of advanced fit ([9]) for a 
spheroid r — 5, k = 5. (Blue) dots - numerical data, solid 
(red) curve - fit by Eq. fit dSJ with B 5.3 x 10"^ Pei « 5.2 
and Pe2 ~ 30.3. Dashed (red) line is BPe^. Here Nu is the 
time-averaged Nu(t) over the period of rotation (determined 
in the simulations), and Nuo is the time-averaged Nu(t) for 
smallest available Pe, e.g. Pe = 0.01. 



IV. ANALYTICAL MODEL OF HEAT 
TRANSFER IN A LAMINAR SHEAR FLOW 

In this section we develop an analytical model for 
Nu(fc, r, p). 



A. Orientational statistics of elongated 
nano-particles 

1. Fokker-Plank equations for orientational PDF 

In order to study the orientational statistics of elon- 
gated nano-particles we consider the Probability Distri- 
bution Function (PDF), V{9, </>, <), which is the probabil- 
ity P{9,(j),t) of finding any particular spheroid with its 
axis of revolution in the interval [9,9 + dO] x [(pjcfi + del)] 
on the unit sphere. The PDF is then defined by 



P{9,(l),t)=V{9,4),t) sin 



d9 



(12) 



It was shown by Burgers [19[ that V{9,(f>,t) satisfies a 
generalized Fokker-Planck equation in the presence of a 
shear: 



dV 
'dt 



= -V • (wV) + DrV^r, 



(13) 



where w = (0,9,4>sm9) is the relative velocity on the 
unit sphere of the axis of revolution for a particle with 
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instantaneous orientation {9, (f>) ignoring all Brownian ef- 
fects. The explicit form of this equation in spherical co- 
ordinates is: 



dr{e, 



1 



dt 



sin 9 



d d ■ 

-F,sin^+-f, 



= 



(14a) 



Here Fq and F^, the 9 and (f> projections of the probability 
density fluxes that have two shear-induced contributions, 
one proportional to w and the other to the rotational 
diffusion coefficient, D^: 



F„ 



89 a 

di^ "89 



F^ = [sine 



V{9,<P) 

90 D, d 
dt s\n9 dcf) 



V{9, 



(14b) 
(14c) 



Burgers in Ref. |l9| derived the rotational diffusion coeffi- 
cient, Dr, of elongated rigid spheroids of revolution: 



1 



r^K3 + Ki 



(15a) 



Here fi is the dynamical fluid viscosity, = 1.374 x 
10~^^erg/deg K is the Boltsmann constant, particle vol- 
ume V = ^TTo^b, 2a = d, b ^ ra, r > 1 and 



rdx 



(r2+2;)3/2(l-Hx) 



(15b) 



- 1 



r arccosh(r) 
- 1 

rdx 



ln(2r/e) 



{r- 
r 



arccosh(r) 



(15c) 



1 



This asymptotics is formally valid for r ^ 1, but give 
better than 10% accuracy already for r > 4, allowing us 
to suggest the approximate "practical" formula 



k^T ln(4rVe) 
4^(r2 +2r/3- 1) ' 



(16) 



that works with an accuracy of about 10% for any r > 1 
and better than with 5%-accuracy for r > 10. 

The complete analysis of Eqs. ([H)) is very involved, 
see e.g. Refs. [l^ and [ll|. We consider only the small- 
diffusion limit, see the next subsection. 



2. Statistics of elongated nano-particles in the small 
diffusion limit 



Jeffery (9| has shown that if inertial and Brownian mo- 
tion affects are completely neglected, then the motion of 
the axis of revolution of a spheroidal particle is described 

by 



tan = r tan r , 



tan9 = CV cos^ r + sin^ r 
= Cr (r^ cos^ 



(17a) 
(17b) 



sin^ 6]-'^\ 



where r = t/T,- with Tj- = {r + l/r)/S, and the constant 
of integration C is called the (Jeffery) orbit constant. 

To analyze the small-diffusion limit we introduce two 
time-scales. The first one defines the periodic motion 
that a nano-particle with a finite r exhibits: Tr = T'j/27r, 
where Tj = 2TT{r + l/r)/S is the Jeffery's period. The 
second time scale is determined by the inverse shear, 
Ts — S~^. Clearly, for large r this is a much shorter 
time scale, so for r ^ 1, 



Tr ~ r/S ~ rTi 



s ■ 



(18) 



It means that the particles spend most of their time near 
pre- and post-aligned states. If the Brownian diffusion 
is small enough such that T^n = 3> Tr one can 

neglect the effect of the Brownian motion on the dynamic 
motions of particles along Jeffery orbits @ even during 
their slow time evolution. In this case the stationary 
Fokker-Plank Eq. takes the simple form 

V • (wV) = . (19) 

This equation can be solved [13, [Hi , giving 

r{9,^)^fiC,r)rc{0,<P), (20a) 

Vc{0,(p) = (r sm9cos^9^cos^(f> + sin^<p/r^ ^ , (20b) 

where Vc{0, 4>) is the PDF along a particular Jeffery or- 
bit with a given integration constant C, and /(C, r) is 
the probability to occupy this orbit. This function is 
normalized as follows: 

V(C)dC = i-. 
/(C, r) was found in [l^ for the limiting cases: 
/(C,r = l) 



/(C,r»l)~- 



/(C,r«l)~- 



1 


c 


47r (C2 


+ 1)3/2 ' 


1 


G 


TT (4C2 


+ 1)3/2 ' 


1 


Cr2 



TT (4r2C2 -I- 1)3/2 



(21a) 
(21b) 
(21c) 



Note that Eq. ()21a|) is exact, providing a consistency 
check of the present approach by comparison with spheri- 
cal particles. The asymptotic, Eq. (j21bp . is very accurate 
for r > 20, but already for r > 4 it provides reasonable 
accuracy (better then ~ 10%). 

The analysis of Eqs. (|2ip together with the available 
numerical solutions for r = 0.26 and 3.86 [l3| allowed us 
to suggest the following approximation 



f{C,r) 



C/(r+\/rY 



3/2 



(22) 



A comparison of this approximation with the exact nu- 
merical solution provided in Ref. [loj is presented in 
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FIG. 7. (Color online). Upper panel: Comparison of "ex- 
act numerical", solutions, obtained from Eq. (|19|l (solid lines) 
and approximate solutions Eq. p2|) for different r (dashed 
lines). There is visible difference (about 5% in the value of 
maximum) only for r = 5. Lower panel: Comparison of 
dependence of (cos^ 6) vs. r obtained by a) applying "exact 
numerical" PDF (black solid line), b) approximate analytical 
PDF, Eq. (|22|l . (red dashed line) and c) analytical approxi- 
mation for this dependence Eq. (|24p (blue dot-dashed line). 

Fig. [71 upper panel. One sees that Eq. ([^ fits the nu- 
merical data with an accuracy better than 5%. This is 
more then enough for our purpose to offer an approxi- 
mate formula for (cos^ 6'n) as a function of r, see below. 

Next we use the fact that JcfFcry orbits do not inter- 
sect on the unit sphere. In other words, fixing C results 
in a relation between 9 and along an orbit. This rela- 
tionship is obtained by inverting Eq. (|17bp : 

C = C{e,(t>)= tan e^r-'^ s\t? + cos2 cf) ■ (23) 

Substituting this function into any one of Eqs. (|2ip wc get 
/(C, r) for a given regime of r. This is then substituted in 
Eqs. (|20|) . leading finally to a solution of the oricntational 
PDF 1^(9, (f>), which can be used for averaging Eqs. pT|) 
in the case of weak Brownian motions. As a consistency 
check of the approach one can consider the trivial case 
r ~ 1. From Eq. (|23p one finds C = tan6', then from 
Rq. (PTa|) /(C) = sin6lcos2 6l/47r and, finally from ((20b)) 



Vc = l/sin0cos^0. As the result one has for sphere 
V = l/47r, as expected. 

For moderate and strong Brownian rotational motion 
the notion of separate JefFery orbits becomes irrelevant. 
In this case we need to solve Eq. (|13|) without approx- 
imations. Once this equation is solved we can com- 
pute (cos^ 9) and substitute the answer in Eqs. (|TT|) . To 
achieve this in the most general case is not a simple task, 
and here we satisfy ourselves with the two limiting cases 
of very large and very small rotational diffusion. The sec- 
ond case was discussed above. For the case of very strong 
rotational diffusion we can use the same Eqs. (fTTj). but 
averaged with a uniform PDF 'P{9,<j)) = I/Att. This is 
because the very strong rotational diffusion tends to dis- 
tribute particles motions around the Jcffcry orbits uni- 
formly. The corrections up to ©(S'/Dj.)^ to such a uni- 
form distribution may be found in Ref. [2^ [2l| 



B. Predictions of the model 

1. Spheroidal particle in a strong Brownian diffusion limit 

With very strong Brownian diffusion, the particles are 
oriented completely randomly, and (cos^ 6'n) = 1/3. In 
this case Eqs. (jlOp and (|TT|) give the results reported in 
Fig. ISl upper panels. The upper left panel shows Nu vs. 
r dependence for various k from 1 to 8600 and oo with 
volume fraction ip = 1%. These results are rather obvi- 
ous: for k ~ 1 one obtains Nu= 1, i.e. no enhancement; 
the larger the k, the larger the heat flux enhancement; 
for any finite k there is a saturation of Nu(r) for ?' — > oo. 
The value of (Nu— 1) may be huge for spheroids (essen- 
tially, rod- like particles at r 3> 1), e.g. for k = 8600 
(diamond in water), Nu— 1 ~ 30, while for spherical par- 
ticles (r = 1), Nu— 1 = 0.03 at the same volume frac- 
tion if ~ 0.01. The values of Nu(fc,r) are bounded, i.e. 
Nu(l,r) < Nu{k,r) < Nu(oo,r). 

The upper right panel shows Nu vs. k dependence at 
(fi = 1% for three values of r: 1, 5 and 10. The values 
of Nu(/c, r) are bounded, too: Nu(fc, 1) < Nu(fc, r) < 
Nu(r, oo). Here, the more elongated particle is (larger 
r), the better the enhancement. 

And last, but not least: elongated particles may touch 
each other much easier. As we show in the Appendix iBl 
the basic geometrical requirement that the mean inter- 
particle distance is less than the largest particle size 
means that r < tq = which is the basic crite- 

ria of the dilute limit. The Nu(r) dependence for r < rp 
is shown by solid curves, while the region r > tq is shown 
by dashed curves in Figs.jSJ left panels. Provided, r < 10 
at ip = 1%, the enhancement may be still considered as 
large but not huge: for fc = 5, the saturation level of 
(Nu-1) is about 0.25 while for spheres Nu-1 sa 0.017. 

The conclusion in the case of very strong Brownian 
diffusion is that the particles with larger k and r bring 
larger heat flux enhancement in the laminar shear flow. 
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2. Spheroidal particle in a weak Brownian diffusion limit 

To complete the calculations of Nu(r) dependence in 
a weak Brownian diffusion limit in the framework of 
model pT|) . we have to find (cos^ 6'n) vs. r. To make 
a long story short, we compared in Fig. [3 lower panel, 
the dependence of (cos^ 0n) vs. r obtained by numerically 
solved PDF, Eq. ^ (see Refs. for more details), 

shown by (black) solid curve, and by approximate ana- 
lytical PDF, Eq. (|22p . shown by (red) dashed curve. As 
one sees these two dependence coincide within line width. 
Moreover, by careful analysis of various limiting cases, we 
suggested the following simple model dependence 



shown in Fig. [71 lower panel, as a (blue) dash-dotted 
curve. Eq. (|24p fits the exact dependence with an ac- 
curacy of 5%. Therefore, Eqs. (|11|) and (|24)) can be 
used in our analysis to make predictions on the ther- 
mal properties in the limit of small Brownian diffusion 
of fluids laden with spheroidal nano-particles of different 
aspect ratios and different thermal conductivity ratios 
with Peclet number up to unity. Corresponding results 
are shown in Fig. [51 lower panels. 

The lower left panel in Fig.[8]shows Nu(r) for different 
k and ip ~ 0.01, at which the Maxwell-Garnctt limit of 
Nu— 1 for spherical particles (r = 1) is 3(/3 = 0.03. Notice, 
the Nu(7') dependence is not monotonic and has a maxi- 
mum at some r — ropt, which depends on k. The reason 
for this is the competition of two effects: more elongated 
nano-particles give larger contribution to the heat flux 
when their longer axis is aligned with the temperature 
gradient, which is orthogonal to the velocity gradient 
(shear) in our case. However longer nano-particles are 
affected more readily by the shear, which tends to ori- 
ent them in the unfavorable direction orthogonal to the 
temperature gradient; then their contribution to the heat 
flux is even less than the one of the spherical particles (cf. 
Fig. [5]). Again, the values of Nu(A:,r) are bounded, i.e. 
Nu(l,r) < Nu(A:,r) < Nu(cxj,r). 

The lower right panel in Fig. [5] shows Nu(fc) for dif- 
ferent aspect ratios, r ai ^ 1%. This is again a 
consequence of the above described competition. The 
values of Nu are bounded by Nu(fc,oo) < Nu(fc,r) < 
Nu[fc, ropt(fc)]. Moreover, for 1 < fc < 7 the optimal 
nano-particle shape is spherical. 

As seen in Fig. [51 lower left panel, there exists a maxi- 
mum of Nu(r ~ fo-pt) for a given k. This maximal Nusselt 
number at its maximizing (optimal) r = ropt(fc) is shown 
in Fig. [9l left, as a function of k ioi ip — 0.01. For k > 7, 
the maximal Nu behaves like 

Nu^*a^ « 1 [0.75 + 0.661n(fc)] , ip = 0.01 .(25a) 

Since here <C 1, the part in parenthesis may be consid- 
ered as (y9-independent, thus, the fit (|25ap may be used 
at any ip < 1%. 



The right panel of Fig.[9]exhibits a contour-density plot 
of Nu(r, fc) at = 0.01. Thick (black) solid and dashed 
curves show ropt(fc) dependence (solid curve is for r < rg, 
and the dashed one is for r > ro = (p'^/^ = 10). The 
inverse dependence fc(ropt) is easy to obtain in a sym- 
bolic computation software by solving 9Nu(r, fc, ip)/dr ~ 
dlC{r,k,ip)/dr = at r = Topt, though, the answer ap- 
pears to be very cumbersome to be shown here. However, 
our analysis reveals that ropt ~ V^, and for ip = 0.01 and 
1 < ropt < 10, . 

r^*t(fc) w -3.03-K 1.57VI, ip = 0.01 . (25b) 

This dependence is shown in Fig. [HI right, as a wide- 
dashed (white) curve, which deviates from the analytical 
solution ropt(fc) for Toot > 10- Again, since <C 1, the 
fit (|25bp may be useco at any (p < 1%. 

In conclusion, we should notice that the rich informa- 
tion about heat flux enhancement, shown in the four pan- 
els of Fig. [51 is just an illustration of the analytical depen- 
dence Nu(iy9, fc, r) given by the analytical expressions (|lip 
and (j24p . This is an important result of our modeling. 

CONCLUSIONS 

We presented a study of the physics of the heat fiux 
in a fiuid laden with nanoparticles of different physical 
properties (shape, thermal conductivity, etc). We de- 
veloped a new analytical model for the effective thermal 
properties of dilute nanofluid suspensions. Our model ac- 
counts for nanoparticle rotation dynamics including the 
fiuid motion around the nanoparticles. We note that our 
model reproduces the classical Maxwell-Garnet model in 
the appropriate static limits. 

We used a combination of theoretical models and nu- 
merical experiments in order to make progress from the 
simplest case of spherical nanoparticles in a quiescent 
fiuid to the most general case of rotating spheroidal par- 
ticles in shear fiows. The new physical ingredient that we 
consider is the exact dynamics of particles in shear flows. 
This constitute a novelty as most of the models intro- 
duced so far, to explain the thermal properties of ther- 
mal colloids, have focused only on the static properties of 
the nanoparticle suspension. Our model starts from the 
realization that particles (spherical or spheroidal) in the 
presence of a gradient of the velocity field are induced to 
rotate. The dynamics of rotation is absolutely non triv- 
ial, but it has been studied at length with correspondece 
(for the case of a laminar and stationary shear flow) to 
the Jeffery orbits. 

The particles rotation dynamics has a double influence 
on the thermal properties of the nanofiuid. First, parti- 
cles rotation induces fluid motions in the proximity of 



^ According to Eqs. mOal l and llll l for ip <^ 1, 9Nu/9r ~ 
d [l3i + {^2 — /3l)(cos■^^?N)] /dr, which is a function of (fc, r) only. 




FIG. 8. (Color online), = 1%, ro = 10. Upper panels: Strong Brownian diffusion limit: (cos^ ^n) — 1/3, lower panels: 
Weak Brownian diffusion limit. Left panels: Nusselt number, Nu, vs. particle aspect ratio, r, for various heat diffusivity 
ratios, k: solid curves for r < rg, dashed curves for r > tq (the region, where the particles may, in principle, start touching 
each other). The shaded (green) area is the region of all possible Nu(fc,r) ai ip — 0.01: the lower bound is given by Nu(l,r), 
the upper one is Nu(oo,r). Right panels: Nu vs. k for various, r. Dotted (black) curve is for r = 1, spherical particles, solid 
(blue) curve - r = 5, dashed (red) curve - r = 10, and dot-dashed (brown) curve - r = 1000, lower panel, or r — >■ oo, upper 
panel. The shaded (green) area is the region of all possible Nu(fc, r) a.t ip — 0.01: the lower bound is Nu[fc, 1] in upper panel and 
Nu(fc, cxj) in lower panel, while the upper bound is Nu(A;,oo) in upper panel and Nu[fc, ropt(fc)] in lower panel, where ropt(fc) 
maximizes Nu at every k, see more in the text. The geometrical constrain imposes 1 < r < max(ropt, 10) for p = 1%. 



the particles, this in turn can enhance the thermal fluxes 
by means of advective motions along the direction of the 
temperature gradient. Second, the JefFery dynamics of 
particles leads to a statistical distribution of particles ori- 
entation that depends on a multitude of parameters, e.g. 
the particle aspect ratio, the shear intensity as well as on 
the intensity of thermal fluctuations. The statistical dis- 
tribution of particle orientation has a dramatic influence 
on the heat flux: an elongated particle oriented along the 
temperature gradient increases the thermal flux, while a 
particle with perpendicular orientation reduces it. 

The statistical orientation of particles can thus pro- 
duce a mixed effects with a non-trivial dependence on 
the particle aspect ratio. More Elongated particles can 
enhance the heat flux because of the stronger contribu- 
tion when properly aligned to the temperature gradient 
but, because of shear, more elongated particles are also 
spending more time in the unfavorable direction (i.e. per- 



pendicular to the temperature flux) thus reducing the 
thermal conductivity of the fluid. 

By means of numerical approximations we are able to 
provide closed expressions for the effective conductivity 
of the fluid under several flow regimes and for several 
physical parameters. Our model considerably extends 
classical models for nanofluid heat transfer, like e.g. the 
one of Maxwell- Garnet, and may help to rationalize some 
of the recent experimental findings. In particular, we 
suggest that experiments should consider more carefully 
measurements performed in quiescent and under flowing 
conditions: the particles dynamics may lead to very dif- 
ferent thermal properties in the two cases. 

Finally, the next steps toward a robust predictive mod- 
els for the heat transfer in nanofluids should include the 
effect of surface (Kapitza) resistance and the effect of 
nanoparticle aggregation. Further it would also be ex- 
tremely important to extend the model to the case of 
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FIG. 9. (Color online). Weak Brownian diffusion limit, = 1%, ro = 10. Left: The largest reachable Nusselt number, 
Nu(fc, r-opt), solid (black) curve (ropt maximizes Nu for a given k). Dashed (red) line is the fit by Nu^ax ~ 1 + 1^(0.75 + 0.66 ln(A;)]. 
Insert: ropt vs. k, solid (blue) curve, and the fit r^pj ~ — 3.03 + 1.57\/fc for 1 < ropt < 10, dashed (red) line. Since ip = 0.01 ^ 1, 
this fit may be considered as </3-independent for < 1%. Right: Contour-density plot of Nu(fc, r) color-codded with the 
"temperature" map: the larger the value of Nu, the wormer the color; contours show iso-values of Nu. Thick (black) solid and 
dashed curves represent those pairs of (fc, r) for which Nu is maximal (i.e. solid (blue) curve in the Insert); solid curve is for 
r < ro, and the dashed one is for r > ro. Wide-dash (white) curve is the ropt(fc) fit. 



heat flux in turbulent nanofluids as this case is very rele- 
vant to many applications. In the presence of turbulence 
a particular attention should also be paid to the effect on 
the drag induced by the presence of spherical, rod-like or 
maybe even dcformablc nanoparticlc inclusions. 
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Appendix A: Numerical approach 

The numerical simulation of the conjugated heat trans- 
fer problem, equations (|4all4cp . is performed by means of 
two coupled D3Q19 Lattice Boltzmann (LB) equations 
under the BGK approximation [l3| (for velocity and tem- 
perature fields) and Molecular-Dynamics simulations (for 
particles motion): 



fiix + ci,t + l) - fiix,t) 



gi{x + ci,t + l) - gi{x,t) = 



fi^x^t)- fi{x,t) 



gn^'t)- gi{x,t) 



(Ala) 
, (Alb) 



where fi{x, t) is the Lattice Boltzmann distribution func- 
tion for particles at {x^t) with velocity q (with I = 
0, . . . , 18 for D3Q19), and fi'^ is its equilibrium distribu- 
tion; gi{x,t) is the distribution functions associated with 
the temperature and g'^'^ is its equilibrium distribution. 



The first LB, Eq. (jAlap . evolves the fluid flow outside 
of the rigid particle and its momentum is coupled with 
the particle boundaries by means of a standard scheme, 
as proposed by Ladd [23. The second LB, Eq. (jAlb|) . 



evolves the temperature field, treated as a passive scalar 
as proposed in [23|, solving the conjugated heat trans- 



fer problem simply by means of adjusting the thermal 
conductivity to the correct values in the fluid and inside 
the particle [eqn. (Pa| and (|4b| ]. Thermal and velocity 
boundary conditions, at the top and bottom walls, im- 
pose the LB populations to equal the equilibrium popula- 
tions (corresponding to the desired velocity and tempera- 
ture). This approach can produce small temperature and 
velocity slip which are kept into account by measuring the 
effective temperature and velocity profiles, thus increas- 
ing the accuracy. The code employed is fully parallelized 
by means of MPI libraries [24 1 thus allowing large sys- 
tem sizes, important to study the influence of finite size 
effects. 

Density, momentum and temperature are defined lo- 
cally at {x,t) as coarse-grained (in velocity space) fields 
of the distribution functions 



= ^ /, , pfUt = J2 (^ifi - PtTf = ^ .9; • (A2) 



A Chapman-Enskog expansio n l25i around the local equi- 
libria fi'^{u,p) and g'p-{u,T) [2g leads to the equations 
for temperature and momentum (j4a|) - (|4b[) : the stream- 
ing step on the left hand side of (|Ala[) reproduces the 
inertial terms in the hydrodynamical equations, whereas 
the diffusive terms (dissipation and thermal diffusion) are 
closely connected to the relaxation (towards equilibrium) 
properties in the right hand side, with v and x related to 
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FIG. 10. Angular velocity vs. time. Solid (red) curve - 
theory of Jeffery Eqs. HA4p . for an infinite domain with 
a constant simple shear at infinity and the creeping (Re— 0) 
flow around the particle, dots - numerical results for small 
but finite Re= 0.05. Configuration: 64^, Pr = 10, fc = 1, 
r = 2 (a = 14, & = 7), Pe = 0.5. 



The angular velocity and the period of such a motion 
were predicted by Jeffery for a case of a creeping flow 
around the particle, i.e. when Re—?' 0: 



Sr/{r + l/r) 



u^. = " " ^ , (A4a) 

cos2(27rt/rj) + r2sin2(27rt/rj) ^ ' 

Tj = 27r(r + l/r)/5 . (A4b) 
For non-vanishing Re, the actual period of rotation of 
the particle deviates from that predicted by Jeffery. In 
Fig. [TU] we compare the results of our LBM simulations 
with Eq. (|A4a|) for the case of r = 2 and Re= 0.05. 
There is a 6% difference in periods due to a finite Re, 
but also due to influence of the other particles present 
due to the periodic b.c, and also the influence of the 
walls {2a/ H ~ 0.44). 



the relaxation times rf, Tg [13[. 



Appendix B: Approximation of non-interacting 
particles 



a. Conjugated heat transfer for a spherical particle at rest 

Consider a spherical particle of radius R immersed in a 
quiescent fluid, in which a constant temperature gradient, 
G, is maintained. The temperature distribution is [l^ : 



The approximation of non-interacting spheroids is 
comming from simple geometrical considerations, and it 
is valid provided the particle aspect ratio r < (p~^/^, as 
conflrmed by the following derivation: 



-G p, 



l-k (R 
2 + k \ p 



r<R, (A3) 
Gp, r>R, 



where p = x'^ + y'^ + z"^ is the distance from the parti- 
cle's center. 

In our Couette flow simulations, G = y A/i?, but the 
temperature boundary conditions are different from [TJ] ■ 
One should expect then deviations close to the walls es- 
pecially for large particles. The results are presented in 

Fig.m 



b. Spheroid in a shear flow - "Jeffery" test 

A spheroidal particle of aspect ratio r in a simple shear 
flow undergoes a spinning motion (in the XY -plane). 



^^11 pariclcs ^ ^^articlo ^^articlc 



Vtotal Vtotal VcS 



Vr, 



off 



£3 



£3 



C = \ = a 



2a < C 



r < 



r < Lp' 



-1/2 



(Bl) 



Here N is the number of particles in the total volume (of 
fluid and particles together), Vtotai, V^ff = Vtotai/A^ is an 
"effective" volume per particle, and C = V^^^ is a char- 
acteristic length/dimension of the effective box/volume 
embedding the particle. 
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